library(dplyr)
library(tibble)
library(DBI)
library(sf)
library(gimms)
library(ncdf4)
library(rgdal)
library(stringr)
library(ggplot2)
library(RColorBrewer)
library(colorspace)
library(cowplot)
library(png)
library(spdep)
library(spatialreg)
library(kableExtra)
library(reshape2)
library(grid)
library(gridExtra)
library(ggpubr)
source("f_plotspatial.R")
memory.limit(size = 160000)
## [1] 160000
library(captioner)
fig_nums <- captioner()
table_nums <- captioner(prefix = "Table")
The coastal drainage basins polygons were computed with a seamless digital elevation model for Fennoscandia (Finstad 2017). The model resulted in more than 200 000 coastal drainage basins, most of them being extremely small polygons covering the coastal regions. In order to simplify the model and reduce the amount of data, we keep only the 0,001 % of the coastal drainage basins that are the biggest. They account to more than 93% of the total surface covered by the initial amount of polygons.
con <- dbConnect(RPostgreSQL::PostgreSQL(),user = "camille.crapart", password = "camille",host = "vm-srv-wallace.vm.ntnu.no", dbname = "nofa")
query <- paste("SELECT gid,geom FROM ","\"Hydrography\"",".waterregions_dem_10m_nosefi",sep = "")
query.2 <- paste("SELECT * FROM ","\"Hydrography\"",".Fenoscandia_LakeRegister",sep = "")
water.regions <- st_read(con,query = query)
lake.register <- st_read(con,query = query.2)
saveRDS(water.regions,"water.regions.Rdata")
water.regions <- readRDS("water.regions.Rdata")
ggplot()+geom_sf(data=water.regions,aes(fill = gid))
water.regions <- readRDS("water.regions.Rdata")
water.regions$area <- st_area(water.regions) %>% as.numeric()
area1 <- sum(water.regions$area)
wr.quantile <- quantile(water.regions$area, probs = c(0.99,0.995,0.999))
wr.sf <- water.regions %>% filter(area > wr.quantile[3])
area2 <- sum(wr.sf$area)
ratio <- area2/area1*100
map.wr <- ggplot(wr.sf) + geom_sf(aes(fill=gid,col=gid))
ggsave(plot = map.wr, filename = "WR/wr.sf.png",dpi = "retina", width = 10, height = 15, units = "cm")
saveRDS(wr.sf,"WR/wr.sf.rds")
st_write(wr.sf, "wr.sf.shp")
NDVI values are extracted from the GIMMS NDVI3g dataset (The National Center for Atmospheric Research 2018), stored on http://poles.tpdc.ac.cn/en/data/9775f2b4-7370-4e5e-a537-3482c9a83d88/, in the “ecocast” dataset and accessed via the “gimms” package (Detsch 2021) on Rstudio.
Data is taken bi-monthly (Tucker et al. 2005). Raster layers of all slices are downloaded using the gimmsdownload function, then monthly composite are calculated using the monthlyComposites function. The maximum NDVI value is taken for each pixel. Afterwards the NDVI values for each catchment polygon are extracted using the extract function from the raster package (Hijmans et al. 2018).
The values stored in the ecocast dataset are composed of the NDVI value and a flag value indicating the goodness of the data (Detsch 2021). All the NDVI values extracted, corresponding to the values for the studied catchments, had a flag of 1, indicating a good value. The NDVI value was retrieved from the stored value using the formula \(NDVI = floor(ndvi3g/10)/1000\) (Detsch 2021), and then the mean of the 3 summer values (June, July and August) was computed for each catchment.
ndvi.max <- readRDS("ndvi.max.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
summer.mean <- raster::stack(ndvi.max[[6]],ndvi.max[[7]],ndvi.max[[8]]) %>% mean()
summer.ndvi <- reclassify(summer.mean,cbind(-Inf, 0, NA), right=FALSE)
summer.scandinavia <- raster::crop(summer.ndvi,c(0,35,55,73))
writeRaster(summer.mean,"ndvi1994.tif",overwrite = T)
ndvi.wr <- raster::extract(summer.ndvi,wr.sf,fun = mean, df = T, sp = T, na.rm = T)
ndvi.wr$area <- NA
ndvi.wr$ndvi <- (floor(ndvi.wr$layer)/10)/1000
saveRDS(ndvi.wr, "WR/ndvi.wr.rds")
knitr::include_graphics("WR/wr.map.NDVI.png")
The runoff data was downloaded from the CORDEX database: https://portal.enes.org/data/enes-model-data/cordex. The data is available at https://esg-dn1.nsc.liu.se/projects/esgf-liu/.
The names of the variables in the CORDEX project follow the CF Metadata convention: http://cfconventions.org/
THe mean of the 30 years between 1970 and 2000 is used (p11 on EURO-CORDEX guidelines https://www.euro-cordex.net/imperia/md/content/csc/cordex/guidance_for_euro-cordex_climate_projections_data_use__2021-02_1_.pdf). This time period matches the temperature and precipitation data provided by WorldClim.
runoff.raster <- "CORDEX/mrros_Lmon_CNRM-CM6-1-HR_historical_r1i1p1f2_gr_185001-201412.nc"
runoff.nc <- nc_open(runoff.raster)
runoff.stack <- raster::stack(runoff.raster, bands = c(1441:1812))
runoff.mean <- runoff.stack %>% mean()
wr.sf <- readRDS("WR/wr.sf.rds")
runoff.wr95 <- raster::extract(runoff.mean, wr.sf, fun = mean, df = T, sp = T, na.rm = T)
names(runoff.wr95) <- c("gid","area","runoff")
saveRDS(runoff.wr95,"WR/runoff.wr95.rds")
knitr::include_graphics("WR/wr.map.Runoff.png")
The Land Cover data was downloaded from https://land.copernicus.eu/pan-european/corine-land-cover/lcc-2000-2006. The previous version (1995-2000) excludes Norway so the 2000 version was preferred.
The land cover data was downloaded as a raster (.tiff file), which included the legend recording the code for each land use. The raster and legend were assembled in QGIS before being imported as a shapefile in R. The codes (corresponding to the line in the extracted dataframe) of the categories of interest were:
Arable land:
Bogs:
Forest:
Bare:
wr.sf <- readRDS("WR/wr.sf.rds")
corine <- raster::raster("Corine_Land_Cover/U2006_CLC2000_V2020_20u1.tif")
wr.corine <- st_transform(wr.sf, crs(corine))
extract_corine <- function(x){
clc <- extract(corine,x)
saveRDS(clc,file = paste("WR/corine/clc",x$gid,"rds",sep="."))
}
wr.list <- split(wr.corine,seq(nrow(wr.corine)))
lapply(wr.list[1143:1440],extract_corine)
#clc.test <- readRDS("WR/clc/clc.16711649.rds")
clc.files <- list.files(path = "WR/corine/", pattern = "clc.*.rds", full.names = T)
list.clc <- sapply(clc.files,readRDS)
list.clc.wr <- sapply(list.clc, as.integer)
names.list <- str_extract(clc.files,"\\d+")
names(list.clc.wr) <- names.list
saveRDS(list.clc.wr,"WR/corine/list.clc.wr.rds")
list.clc.wr <- readRDS("WR/corine/list.clc.wr.rds")
clc.tab.wr <- sapply(list.clc.wr, function(x) tabulate(x,45)) %>% t()
clc.tab.area.wr <- clc.tab.wr*prod(res(corine))
catchment.area <- rowSums(clc.tab.area.wr)
clc.tab.prop.wr <- sweep(clc.tab.area.wr,1,catchment.area, FUN = "/")
saveRDS(clc.tab.prop.wr,"WR/corine/clc.tab.prop.wr.rds")
clc.tab.prop.wr <- readRDS("WR/corine/clc.tab.prop.wr.rds") %>% as.data.frame() %>% tibble::rownames_to_column(var = "gid")
bogs <- clc.tab.prop.wr %>% dplyr::select(c("gid","V36")) %>% setNames(c("gid","bogs"))
saveRDS(bogs,"WR/corine/bogs.wr.rds")
arable <- data.frame( gid = clc.tab.prop.wr$gid, arable = rowSums(select(clc.tab.prop.wr, c("V12","V16","V18","V19","V20"))))
saveRDS(arable,"WR/corine/arable.wr.rds")
forest <- data.frame( gid = clc.tab.prop.wr$gid, forest = rowSums(select(clc.tab.prop.wr, c("V23","V24","V25"))))
saveRDS(forest,"WR/corine/forest.wr.rds")
bare <- data.frame( gid = clc.tab.prop.wr$gid, bare = rowSums(select(clc.tab.prop.wr, c("V31","V32","V33"))))
saveRDS(bare,"WR/corine/bare.wr.rds")
water <- data.frame( gid = clc.tab.prop.wr$gid, water = rowSums(select(clc.tab.prop.wr, c("V40","V41","V42","V43","V44"))))
saveRDS(water,"WR/corine/water.wr.rds")
list.grob <- c("WR/wr.map.Bog.png", "WR/wr.map.Arable.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow = 1, ncol = 2, labels = c("a)","b)"), label_size = 10)
Temperature and precipitation are downloaded from the WorldClim database: https://worldclim.org/data/index.html
Temperature in historical Worldclim is in °C * 10 (to reduce file size). Precipitation are average monthly precipitation in mm. Here we use the mean monthly precipitation and temperature for the period 1970-2000.
wr.sf <- readRDS("WR/wr.sf.rds")
bio.fennoscandia <- raster::getData(name = "worldclim", var = "bio", res = 2.5)
temp.wr <- raster::extract(bio.fennoscandia$bio1, wr.sf, fun = mean, df = T, sp = T, na.rm = T)
names(temp.wr) <- c("gid","area","temp")
temp.wr$temp <- temp.wr$temp / 10
saveRDS(temp.wr,"WR/temp.wr.rds")
#temp.wr <- readRDS("WR/temp.wr.rds")
prec.wr <- raster::extract(bio.fennoscandia$bio12, wr.sf, fun = mean, df = T, sp = T, na.rm = T)
names(prec.wr) <- c("gid","area","prec")
saveRDS(prec.wr,"WR/prec.wr.rds")
#prec.wr <- readRDS("WR/prec.wr.rds")
tp <- cbind(wr.sf[,1],temp.wr[,3]) %>% cbind(prec.wr[,3])
saveRDS(tp,"WR/tp.rds")
list.grob <- c("WR/wr.map.Temp.png", "WR/wr.map.Precip.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow= 1, ncol = 2, labels = c("a)", "b)"), label_size = 10)
The nitrogen deposition data was extracted from the EMEP database (https://emep.int/mscw/mscw_moddata.html#Comp). The data from 2000 was used as it is the first model available.
N-deposition (NOx and NH3) is the sum of:
library(ncdf4)
EMEP_file <- "EMEP/EMEP01_rv4.42_year.2000met_2000emis_rep2021.nc"
wr.sf <- readRDS("WR/wr.sf.rds")
woxn <- raster(EMEP_file, varname = "WDEP_OXN")
doxn <- raster(EMEP_file, varname = "DDEP_OXN_m2Grid")
wrdn <- raster(EMEP_file, varname = "WDEP_RDN")
drdn <- raster(EMEP_file, varname = "DDEP_RDN_m2Grid")
woxn_wr <- extract(woxn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(woxn_wr) <- c("gid","area","woxn")
saveRDS(woxn_wr,"woxn_wr.rds")
doxn_wr <- extract(doxn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(doxn_wr) <- c("gid","area","doxn")
saveRDS(doxn_wr,"doxn_wr.rds")
wrdn_wr <- extract(wrdn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(wrdn_wr) <- c("gid","area","wrdn")
saveRDS(wrdn_wr,"wrdn_wr.rds")
drdn_wr <- extract(drdn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(drdn_wr) <- c("gid","area","drdn")
saveRDS(drdn_wr,"drdn_wr.rds")
ndep.wr <- merge(woxn_wr,doxn_wr, by = "gid") %>% merge(wrdn_wr, by = "gid") %>% merge(drdn_wr, by = "gid")
ndep.wr$tndep <- rowSums(ndep.wr@data[,c(3:6)])
saveRDS(ndep.wr,"WR/ndep.wr.rds")
knitr::include_graphics("WR/wr.map.TNdep.png")
The altitude data was accessed through the “getData” function, and comes from the NASA digital elevation model (SRTM 90 m). The slope was then computed using the “terrain” function from the raster package.
wr.sf <- readRDS("WR/wr.sf.rds")
alt.nor <- raster::getData("alt",country = "NOR", mask = T)
alt.swe <- raster::getData("alt",country = "SWE", mask = T)
alt.fin <- raster::getData("alt",country = "FIN", mask = T)
alt.list <- list(alt.nor,alt.swe,alt.fin)
alt <- do.call(raster::merge,alt.list)
alt.wr <- raster::extract(alt, dplyr::select(wr.sf, c("gid","geom")), sp = T, fun = mean, df = T, na.rm = T)
names(alt.wr) <- c("gid","alt")
saveRDS(alt.wr,"WR/alt.wr.rds")
slope.nor <- terrain(alt.nor, unit = "degrees")
slope.swe <- terrain(alt.swe, unit = "degrees")
slope.fin <- terrain(alt.fin, unit = "degrees")
slope.list <- list(slope.nor,slope.swe,slope.fin)
slope <- do.call(raster::merge,slope.list)
slope.nona <- reclassify(slope,cbind(NA,100), right = NA)
slope.wr <- raster::extract(slope,wr.sf, sp = T, df = T, fun = mean, na.rm = T)
names(slope.wr) <- c("gid","area","slope")
saveRDS(slope.wr, "WR/slope.wr.rds")
The data was gathered and merged into a single dataframe, based on the catchment identifyer (“gid”).
ndvi.wr <- readRDS("WR/ndvi.wr.rds")
ndvi.wr$area <- NULL
runoff.wr95 <- readRDS("WR/runoff.wr95.rds")
runoff.wr95$area <- NULL
bogs.wr <- readRDS("WR/corine/bogs.wr.rds")
forest.wr <- readRDS("WR/corine/forest.wr.rds")
arable.wr <- readRDS("WR/corine/arable.wr.rds")
bare.wr <- readRDS("WR/corine/bare.wr.rds")
water.wr <- readRDS("WR/corine/water.wr.rds")
#clc.wr.sum <- readRDS("WR/clc/clc.wr.sum.Rdata")
alt.wr <- readRDS("WR/alt.wr.rds")
slope.wr <- readRDS("WR/slope.wr.rds")
slope.wr$area <- NULL
tp <- readRDS("WR/tp.rds")
ndep.wr <- readRDS("WR/ndep.wr.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
all.wr.sf <- merge(wr.sf,runoff.wr95,by = "gid") %>% merge(ndvi.wr,by="gid") %>% merge(bogs.wr,by = "gid") %>% merge(forest.wr, by = "gid") %>% merge(arable.wr, by = "gid") %>% merge(bare.wr, by = "gid") %>% merge(water.wr, by = "gid") %>% merge(st_drop_geometry(tp),by = "gid") %>% merge(alt.wr, by = "gid") %>% merge(slope.wr, by = "gid") %>% merge(ndep.wr, by = "gid")
saveRDS(all.wr.sf, "WR/all.wr.sf.rds")
The data was then cleaned and the centroid of each coastal drainage basin was computed.
all.wr.sf <- readRDS("WR/all.wr.sf.rds")
wr.sf.95 <- all.wr.sf %>% dplyr::filter(ndvi != 0) %>% filter(runoff > 0) %>% filter(is.na(bogs) == F) %>% filter(is.na(alt) == F) %>% filter(is.na(slope) == F)
wr.sf.95 <- wr.sf.95[!duplicated(wr.sf.95$gid),]
wr.sf.95$Runoff <- wr.sf.95$runoff *365*24*60*60 # from kg/m2/s to kg/m2/year or mm/y
wr.sf.95$logRunoff <- log(wr.sf.95$Runoff)
sfc_point_to_matrix = function(x) {
matrix(
unlist(x, use.names = FALSE),
nrow = length(x),
byrow = TRUE,
dimnames = list(1:length(x))
)
}
wr.centroid <- st_centroid(wr.sf.95, of.largest.polygon = T)
coord <- sfc_point_to_matrix(wr.centroid$geometry) %>% as.data.frame() %>% setNames(c("Longitude","Latitude"))
wr.sf.95 <- cbind(wr.sf.95, coord)
wr.sf.95$area <- st_area(wr.sf.95)
saveRDS(wr.sf.95, "WR/wr.sf.95.rds")
noforest <- filter(wr.sf.95,forest == 0) #%>% dplyr::select(noforest,c("gid","area"))
#ggplot(noforest[which.max(noforest$area),])+geom_sf(aes(fill=as.factor(gid)))
weird.gid <- noforest$gid[which.max(noforest$area)]
weird.cat <- wr.sf.95 %>% filter(gid == weird.gid) %>% st_transform(CRS("EPSG:3035"))
corine <- raster::raster("Corine_Land_Cover/U2006_CLC2000_V2020_20u1.tif")
#weird.clc <- raster::extract(corine,weird.cat)
#saveRDS(weird.clc,"WR/clc/weird.clc.rds")
weird.clc <- readRDS("WR/clc/weird.clc.rds")
weird.tab <- sapply(weird.clc, function(x) tabulate(x,45))
weird.tabl.area <- weird.tab*prod(res(corine))
catchment.area <- colSums(weird.tabl.area)
weird.tab.prop <- sweep(weird.tabl.area,2,catchment.area, FUN = "/")
wr.sf.95$bogs[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[36]
wr.sf.95$forest[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[24]+weird.tab.prop[24]+weird.tab.prop[25]
wr.sf.95$arable[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[12]+weird.tab.prop[16]+weird.tab.prop[18]+weird.tab.prop[19]+weird.tab.prop[20]
wr.sf.95$glacier[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[34]
wr.sf.95$bare[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[31]+weird.tab.prop[32]+weird.tab.prop[33]
saveRDS(wr.sf.95,"WR/wr.sf.95.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
toreplace <- which(names(wr.sf.95) %in% c("ndvi","bogs","forest","arable","bare","water","temp","prec","alt","slope","tndep"))
names(wr.sf.95)[toreplace] <- c("NDVI","Bog","Forest","Arable","Bare","Water","Temp","Precip","Alt","Slope","TNdep")
wr.sf.95$area.x <- NULL
wr.sf.95$area.y <- NULL
saveRDS(wr.sf.95, "WR/wr.sf.95.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
m <- c("NDVI","Runoff","Precip","Temp","Bog","Forest","Arable","TNdep")
u <- c("NDVI on scale 0 to 1", "Mean Runoff, mm/y","Mean monthly precipitation,\n mm","Mean annual temperature,\n °C","Proportion of bogs","Proportion of Forest","Proportion of\narable land","Total Nitrogen Deposition (mg/m2)")
couleurs <- c("Greens","Blues","Blues","RdYlBu","YlOrBr","RdYlGn","RdYlGn","PuRd")
dir <- c(T,F,F,T,T,F,T,F)
for(i in m){
mysf <- wr.sf.95
map <- ggplot(mysf)+geom_sf(aes(fill = .data[[i]], col = .data[[i]]))+
scale_fill_continuous_divergingx(palette = "RdYlBu",
aesthetics = c("colour","fill"),
rev = dir[grep(i,m)],
name = paste(u[grep(i,m)],"\n(",i,")",sep=""))+
theme_minimal(base_size = 9)+
theme(legend.position = "bottom")
filename_i <- paste("WR/wr.map",i,"png",sep = ".")
ggsave(plot = map, filename = filename_i, dpi = "retina", width = 10, height = 12, units = "cm")
}
wr.sf.95 <- readRDS("WR/wr.sf.95.rds") %>% st_drop_geometry()
wr.sf.95$area_km2 <- wr.sf.95$area * 10^(-6)
wr.min <- lapply(wr.sf.95, min) %>% as.data.frame()
wr.max <- lapply(wr.sf.95, max) %>% as.data.frame()
wr.median <- lapply(wr.sf.95, median) %>% as.data.frame()
wr.mean <- lapply(wr.sf.95, mean) %>% as.data.frame()
wr.sd <- lapply(wr.sf.95, sd) %>% as.data.frame()
wr.sum <- rbind(wr.min, wr.max, wr.median, wr.mean, wr.sd)
rownames(wr.sum) <- c("Min","Max","Median","Mean","Sd")
wr.sum %>% dplyr::select(c("NDVI", "Forest", "Arable","Temp","Precip","TNdep","Runoff","area_km2")) %>% kable(digits = 2) %>%
kable_styling(bootstrap_options = "bordered")
| NDVI | Forest | Arable | Temp | Precip | TNdep | Runoff | area_km2 | |
|---|---|---|---|---|---|---|---|---|
| Min | 0.11 | 0.00 | 0.00 | -4.24 | 434.55 | 63.88 | 55.17 | 27.3728 [m^2] |
| Max | 0.90 | 0.96 | 0.92 | 8.13 | 2842.75 | 2387.52 | 971.21 | 49832.2278 [m^2] |
| Median | 0.73 | 0.38 | 0.01 | 3.24 | 671.00 | 432.82 | 250.68 | 76.2258 [m^2] |
| Mean | 0.68 | 0.39 | 0.08 | 3.05 | 966.33 | 586.44 | 339.81 | 751.3046 [m^2] |
| Sd | 0.15 | 0.25 | 0.15 | 2.97 | 582.82 | 491.72 | 222.63 | 3752.6959 [m^2] |
We model the mean TOC concentration in all the water region in 1995 based on the simple SEM model of the previous section, using NDVI, bogs and runoff as predictors.
The neighbor matrix of the coastal drainage basins polygons were computed using their centroids.
wr.sf.95 <- readRDS("wr.sf.95.rds")
wr.centroid <- st_centroid(wr.sf.95, of_largest_polygon = T)
wr.neigh.set <- knearneigh(wr.centroid, k = 100)
wr.neigh.nb <- knn2nb(wr.neigh.set)
wr.kmat <- nb2listw(wr.neigh.nb)
saveRDS(wr.kmat,"WR/wr.kmat.rds")
The TOC concentration in coastal drainage basins was modeled using the Spatial Error Linear Model with 5 predictors:
library(spatialreg)
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred <- predict(sem.fennoscandia.5,wr.sf.95,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
sem.pred.wr <- cbind(wr.sf.95,wr.pred$fit)
sem.pred.wr$TOC95 <- exp(sem.pred.wr$wr.pred.fit)
map.fit <- ggplot()+geom_sf(data = st_as_sf(sem.pred.wr), aes(fill = TOC95, col = TOC95))+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1,
name = "Modelled TOC 1995 (mg/L)")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
ggsave(plot = map.fit, filename = "WR/wr.toc.png",dpi = "retina", width = 6, height = 8, units = "cm")
map.log <- ggplot()+geom_sf(data = sem.pred.wr, aes(fill = wr.pred.fit, col = wr.pred.fit))+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1,
name = "Modelled log(TOC)")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
ggsave(plot = map.log, filename = "WR/wr.log.toc.png",dpi = "retina", width = 6, height = 8, units = "cm")
saveRDS(sem.pred.wr,"WR/sem.pred.wr.rds")
knitr::include_graphics("WR/wr.log.toc.png")
The forecast of the temperature and precipitation in the periods 2041-2060 and 2081-2100, for SSP 1-2.6 and SSP 3-7.0, were extracted based on the Future Worldclim dataset.
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
ssp126.cnrm.temp <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp126_2041-2060.tif", band = 1)
ssp370.cnrm.temp <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp370_2041-2060.tif", band = 1)
ssp126.cnrm.temp.2050.crop <- raster::crop(ssp126.cnrm.temp,c(0,35,55,73))
png("worldclim_future/ssp126.temp.2050.png")
plot(ssp126.cnrm.temp.2050.crop)
dev.off()
ssp370.cnrm.temp.2050.crop <- raster::crop(ssp370.cnrm.temp,c(0,35,55,73))
png("worldclim_future/ssp370.temp.2050.png")
plot(ssp370.cnrm.temp.2050.crop)
dev.off()
ssp126.temp.wr <- raster::extract(ssp126.cnrm.temp,select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm =T)
names(ssp126.temp.wr) <- c("gid", "temp126")
saveRDS(ssp126.temp.wr,"WR/ssp126.temp.wr.rds")
ssp370.temp.wr <- raster::extract(ssp370.cnrm.temp,select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.temp.wr) <- c("gid", "temp370")
saveRDS(ssp370.temp.wr,"WR/ssp370.temp.wr.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
raster.ssp126.temp.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp126_2081-2100.tif", band = 1)
raster.ssp370.temp.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp370_2081-2100.tif", band = 1)
ssp126.cnrm.temp.2100.crop <- raster::crop(raster.ssp126.temp.2100,c(0,35,55,73))
png("worldclim_future/ssp126.temp.2100.png")
plot(ssp126.cnrm.temp.2100.crop)
dev.off()
ssp370.cnrm.temp.2100.crop <- raster::crop(raster.ssp370.temp.2100,c(0,35,55,73))
png("worldclim_future/ssp370.temp.2100.png")
plot(ssp370.cnrm.temp.2100.crop)
dev.off()
ssp126.temp.2100 <- raster::extract(raster.ssp126.temp.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm =T)
names(ssp126.temp.2100) <- c("gid", "temp126.2100")
saveRDS(ssp126.temp.2100,"WR/ssp126.temp.2100.rds")
ssp370.temp.2100 <- raster::extract(raster.ssp370.temp.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.temp.2100) <- c("gid", "temp370.2100")
saveRDS(ssp370.temp.2100,"WR/ssp370.temp.2100.rds")
listgrob <- list("worldclim_future/ssp126.temp.2050.png", "worldclim_future/ssp370.temp.2050.png","worldclim_future/ssp126.temp.2100.png","worldclim_future/ssp370.temp.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol = 2, nrow = 2, labels = c("a) SSP126-2050","b) SSP370-2050","c) SSP126-2100","d) SSP370-2100"), label_size = 10)
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
ssp126.cnrm.prec <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp126_2041-2060.tif", band = 12)
ssp370.cnrm.prec <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp370_2041-2060.tif", band = 12)
ssp126.cnrm.prec.crop <- raster::crop(ssp126.cnrm.prec,c(0,35,55,73))
png("worldclim_future/ssp126.2050.prec.png")
plot(ssp126.cnrm.prec.crop)
dev.off()
ssp370.cnrm.prec.crop <- raster::crop(ssp370.cnrm.prec,c(0,35,55,73))
png("worldclim_future/ssp370.2050.prec.png")
plot(ssp370.cnrm.prec.crop)
dev.off()
ssp126.prec.wr <- raster::extract(ssp126.cnrm.prec,select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp126.prec.wr) <- c("gid", "prec126")
saveRDS(ssp126.prec.wr,"WR/ssp126.prec.wr.rds")
ssp370.prec.wr <- raster::extract(ssp370.cnrm.prec,dplyr::select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.prec.wr) <- c("gid", "prec370")
saveRDS(ssp370.prec.wr,"WR/ssp370.prec.wr.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
raster.ssp126.prec.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp126_2081-2100.tif", band = 12)
raster.ssp370.prec.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp370_2081-2100.tif", band = 12)
ssp126.cnrm.prec.2100.crop <- raster::crop(raster.ssp126.prec.2100,c(0,35,55,73))
png("worldclim_future/ssp126.2100.prec.png")
plot(ssp126.cnrm.prec.2100.crop)
dev.off()
png("worldclim_future/ssp370.2100.prec.png")
ssp370.cnrm.prec.2100.crop <- raster::crop(raster.ssp370.prec.2100,c(0,35,55,73))
plot(ssp370.cnrm.prec.2100.crop)
ssp126.prec.2100 <- raster::extract(raster.ssp126.prec.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp126.prec.2100) <- c("gid", "prec126.2100")
saveRDS(ssp126.prec.2100,"WR/ssp126.prec.2100.rds")
ssp370.prec.2100 <- raster::extract(raster.ssp370.prec.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.prec.2100) <- c("gid", "prec370.2100")
saveRDS(ssp370.prec.2100,"WR/ssp370.prec.2100.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
ssp126.temp.wr <- readRDS("WR/ssp126.temp.wr.rds")
ssp370.temp.wr <- readRDS("WR/ssp370.temp.wr.rds")
ssp126.prec.wr <- readRDS("WR/ssp126.prec.wr.rds")
ssp370.prec.wr <- readRDS("WR/ssp370.prec.wr.rds")
ssp126.temp.2100 <- readRDS("WR/ssp126.temp.2100.rds")
ssp370.temp.2100 <- readRDS("WR/ssp370.temp.2100.rds")
ssp126.prec.2100 <- readRDS("WR/ssp126.prec.2100.rds")
ssp370.prec.2100 <- readRDS("WR/ssp370.prec.2100.rds")
# SSP126
df.ssp126 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff")) %>% merge(ssp126.temp.wr, by = "gid") %>% merge(ssp126.prec.wr, by = "gid") %>% st_drop_geometry() %>% setNames(c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff","Temp","Precip"))
df.ssp126$logPrecip <- df.ssp126$Precip %>% log()
saveRDS(df.ssp126, "WR/df.ssp126.rds")
# SSP 370
df.ssp370 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff")) %>% merge(ssp370.temp.wr, by = "gid") %>% merge(ssp370.prec.wr, by = "gid") %>% st_drop_geometry() %>% setNames(c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff","Temp","Precip"))
df.ssp370$logPrecip <- df.ssp370$Precip %>% log()
saveRDS(df.ssp370, "WR/df.ssp370.rds")
# SSP126 2100
df.ssp126.2100 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff")) %>% merge(ssp126.temp.2100, by = "gid") %>% merge(ssp126.prec.2100, by = "gid") %>% st_drop_geometry() %>% setNames(c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff","Temp","Precip"))
df.ssp126.2100$logPrecip <- df.ssp126.2100$Precip %>% log()
saveRDS(df.ssp126.2100, "WR/df.ssp126.2100.rds")
# SSP 370 2100
df.ssp370.2100 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff")) %>% merge(ssp370.temp.2100, by = "gid") %>% merge(ssp370.prec.2100, by = "gid") %>% st_drop_geometry() %>% setNames(c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff","Temp","Precip"))
df.ssp370.2100$logPrecip <- df.ssp370.2100$Precip %>% log()
saveRDS(df.ssp370.2100, "WR/df.ssp370.2100.rds")
listgrob <- list("worldclim_future/ssp126.2050.prec.png", "worldclim_future/ssp370.2050.prec.png","worldclim_future/ssp126.2100.prec.png","worldclim_future/ssp370.2100.prec.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol = 2, nrow = 2, labels = c("a) SSP126-2050","b) SSP370-2050","c) SSP126-2100","d) SSP370-2100"), label_size = 10)
We model NDVI based on temperature and precipitation, and use this model to predict the future NDVI in 2050. Choice of future data: ssp 126 as best scenario, ssp 370 as worse scenario (Hausfather 2019; CORDEX 2021). CRNM according to Raju and Kumar (2020).
Using a SELM model for NDVI, we notice that high residuals are observed in catchments where original NDVI is lower than 0,4. Myers-Smith et al. (2020) show that the relationship between NDVI and the Arctic biomass is not linear: under 0,4, an increase in NDVI reflects the vegetalization of bare ground, while over 0,4, an increase of NDVI reflects an increased productivity in established forests.
Several models were tested to model the NDVI: a linear model, a spatial error linear model, a polynomial model, and a polynomial model combined with a binomial model. The NDVI was fitted on the 1995 data from the Northern Lakes Survey and tested on on the coastal drainage basins.
library(betareg)
library(mgcv)
library(VGAM)
wr.kmat <- readRDS("WR/wr.kmat.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
fennoscandia <- readRDS("fennoscandia.33.rds")
fennoscandia$Latitude <- fennoscandia$latitude
fen.kmat <- readRDS("k.neigh.w.rds")
ndvi.cor <- cor(dplyr::select(st_drop_geometry(fennoscandia), c("NDVI", "Temp", "Precip", "Arable", "Bare", "Latitude", "Slope","Runoff")))
corrplot::corrplot(ndvi.cor, method = "number")
par(mfrow = c(2,3))
plot(NDVI~Temp+log(Precip)+Arable+Bare+Latitude+Slope, data = st_drop_geometry(fennoscandia))
par(mfrow = c(3,2))
#LM
summary(m1 <- lm(NDVI ~ Temp + log(Precip), data = fennoscandia))
summary(m1 <- step(m1))
plot(predict(m1), fennoscandia$NDVI, pch = ".", main="a ) NDVI - LM", ylab="Observed NDVI", xlab="Predicted",abline(0,1))
text(x = 0.8, y = 0.4, label = paste("r2 =", round(summary(m1)$adj.r.squared,2), sep= ""))
#SELM
fen.kmat <- readRDS("k.neigh.w.rds")
selm.model <- errorsarlm(formula = NDVI~Temp+logPrecip, data = fennoscandia, listw = fen.kmat)
saveRDS(selm.model, "WR/selm.ndvi.rds")
selm.model <- readRDS("WR/selm.ndvi.rds")
AIC <- 2*(3-1)-2*ndvi.model$LL
plot(ndvi.model$fitted.value, fennoscandia$NDVI, pch = ".", main="b) NDVI - SELM", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.8, y = 0.5, label = paste("r2 =", round(cor(ndvi.model$fitted.values, fennoscandia$NDVI)^2,2), sep= ""))
# quadratic
summary(m2 <- lm(NDVI ~ (Temp + I(Temp^2)) * (log(Precip) + I(log(Precip)^2)), data = fennoscandia))
summary(m2 <- step(m2))
plot(predict(m2), fennoscandia$NDVI, pch = ".", main="c) NDVI - quadratic model", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.7, y = 0.4, label = paste("r2 =", round(summary(m2)$adj.r.squared,2), sep= ""))
# poly
summary(m3 <- lm(NDVI ~ poly(cbind(Temp, logPrecip), degree = 2), data = fennoscandia))
summary(m3 <- step(m3))
plot(predict(m3), fennoscandia$NDVI, pch = ".", main="d) NDVI - polynomial ", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.7, y = 0.4, label = paste("r2 =", round(summary(m3)$adj.r.squared,2), sep= ""))
saveRDS(m3,"WR/poly.ndvi.rds")
# glm
summary(m4 <- glm(NDVI ~ poly(cbind(Temp, logPrecip), degree = 2), data = fennoscandia, family = "quasibinomial"))
#summary(m4 <- step(m4))
plot(predict(m4, type = "response"), fennoscandia$NDVI, pch = ".", main="e) NDVI - glm ", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.7, y = 0.4, label = paste("r =", cor(fennoscandia$NDVI, m4$fitted.values), sep= ""))
saveRDS(m4, "WR/glm.ndvi.rds")
#betareg
summary(m5 <- betareg(NDVI ~ poly(cbind(Temp, logPrecip), degree = 2), data = fennoscandia, link = "logit"))
plot(predict(m5, type = "response"), fennoscandia$NDVI, pch = ".", main="e) NDVI - glm ", ylab="Observed", xlab="Predicted",abline(0,1))
saveRDS(m5, "WR/betareg.ndvi.rds")
The SELM and the polynomial model had the best results. We tested them on the coastal drainage basins dataset.
k.neigh.w <- readRDS("k.neigh.w.rds")
fennoscandia <- readRDS("fennoscandia.rds")
fennoscandia.33 <- readRDS("fennoscandia.33.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- log(wr.sf.95$Precip)
selm.ndvi <- readRDS("WR/selm.ndvi.rds")
selm.ndvi.fen <- predict(selm.ndvi,wr.sf.95,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.selm.ndvi.model <- wr.sf.95 %>% st_as_sf() %>% mutate(ndvi.fit = selm.ndvi.fen$fit, ndvi.residuals = wr.sf.95$NDVI - selm.ndvi.fen$fit)
wr.selm.ndvi.model$area <- st_area(wr.selm.ndvi.model) %>% as.numeric()
ggplot(wr.selm.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = Latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(wr.selm.ndvi.model$NDVI,wr.selm.ndvi.model$ndvi.fit)^2,2)), x = 0.25, y = 0.5, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
labs(title = "Observed vs predicted NDVI for SELM")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/wr.ndvi.selm.results.png", dpi = "retina", width = 10, height = 8, units = "cm")
wr.residuals.selm <- ggplot()+
geom_sf(data = wr.selm.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "SELM residuals")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")+
labs(title = "Map of residuals for SELM")
ggsave(plot = wr.residuals.selm, filename = "WR/wr.selm.ndvi.residuals.png",dpi = "retina", width = 12, height = 16, units = "cm")
poly.ndvi <- readRDS("WR/poly.ndvi.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- log(wr.sf.95$Precip)
fennoscandia <- readRDS("fennoscandia.33.rds")
poly.ndvi.wr <- predict(poly.ndvi,wr.sf.95, type = "response")
wr.poly.ndvi.model <- wr.sf.95 %>% mutate(ndvi.fit = poly.ndvi.wr, ndvi.residuals = wr.sf.95$NDVI - poly.ndvi.wr)
ggplot(wr.poly.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = Latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(wr.poly.ndvi.model$NDVI,wr.poly.ndvi.model$ndvi.fit),2)), x = 0.75, y = 0.4, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
labs(title = "Observed vs predicted NDVI for polynomial model")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/wr.ndvi.poly.results.png", dpi = "retina", width = 10, height = 8, units = "cm")
wr.residuals <- ggplot()+
geom_sf(data = wr.poly.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "Polynomial residuals")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")+
labs(title = "Map of residuals for Polynomial model")
ggsave(plot = wr.residuals, filename = "WR/wr.poly.ndvi.residuals.png",dpi = "retina", width = 12, height = 16, units = "cm")
glm.ndvi <- readRDS("WR/glm.ndvi.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- log(wr.sf.95$Precip)
glm.ndvi.wr <- predict(glm.ndvi,wr.sf.95, type = "response")
wr.glm.ndvi.model <- wr.sf.95 %>% mutate(ndvi.fit = glm.ndvi.wr, ndvi.residuals = wr.sf.95$NDVI - glm.ndvi.wr)
ggplot(wr.glm.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = Latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(wr.glm.ndvi.model$NDVI,wr.glm.ndvi.model$ndvi.fit),2)), x = 0.8, y = 0.3, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
labs(title = "Observed vs predicted NDVI for binomial model")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/wr.ndvi.glm.results.png", dpi = "retina", width = 10, height = 8, units = "cm")
wr.residuals <- ggplot()+
geom_sf(data = wr.glm.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "glm residuals")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")+
labs(title = "Map of residuals for Binomial model")
ggsave(plot = wr.residuals, filename = "WR/wr.glm.ndvi.residuals.png",dpi = "retina", width = 12, height = 16, units = "cm")
betareg.ndvi <- readRDS("WR/betareg.ndvi.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- log(wr.sf.95$Precip)
betareg.ndvi.wr <- predict(betareg.ndvi,wr.sf.95, type = "response")
wr.betareg.ndvi.model <- wr.sf.95 %>% mutate(ndvi.fit = betareg.ndvi.wr, ndvi.residuals = wr.sf.95$NDVI - betareg.ndvi.wr)
ggplot(wr.betareg.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = Latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(wr.betareg.ndvi.model$NDVI,wr.betareg.ndvi.model$ndvi.fit),2)), x = 0.8, y = 0.3, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
labs(title = "Observed vs predicted NDVI for binomial model")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/wr.ndvi.betareg.results.png", dpi = "retina", width = 10, height = 8, units = "cm")
wr.residuals <- ggplot()+
geom_sf(data = wr.betareg.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "betareg residuals")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")+
labs(title = "Map of residuals for Binomial model")
ggsave(plot = wr.residuals, filename = "WR/wr.betareg.ndvi.residuals.png",dpi = "retina", width = 12, height = 16, units = "cm")
listgrob1 <- list("WR/wr.ndvi.selm.results.png","WR/wr.ndvi.poly.results.png","WR/wr.ndvi.glm.results.png","WR/wr.ndvi.betareg.results.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob1, ncol = 2, nrow = 2, labels = c("a)","b)","c)","d)"))
listgrob2 <- list("WR/wr.selm.ndvi.residuals.png","WR/wr.poly.ndvi.residuals.png", "WR/wr.glm.ndvi.residuals.png","WR/wr.betareg.ndvi.residuals.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob2, ncol= 2, nrow = 2,labels = c("e)","f)","g)","h)"))
selm.ndvi <- readRDS("WR/selm.ndvi.rds")
poly.ndvi <- readRDS("WR/poly.ndvi.rds")
glm.ndvi <- readRDS("WR/glm.ndvi.rds")
betareg.ndvi <- readRDS("WR/betareg.ndvi.rds")
fennoscandia <- readRDS("fennoscandia.33.rds")
xtemp <- seq(min(fennoscandia$Temp), max(fennoscandia$Temp), 0.1)
yprecip <- seq(min(fennoscandia$logPrecip), max(fennoscandia$logPrecip), 0.1)
data.fit <- expand.grid(Temp = xtemp, logPrecip = yprecip)
mtrx3d <- predict(selm.ndvi, newdata = data.fit, type = "response")
mtrx.melt <- cbind(data.fit, mtrx3d)
names(mtrx.melt) <- c("Temp","logPrecip","NDVI")
cp1 <- ggplot(mtrx.melt,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = fennoscandia, aes(x=Temp, y = logPrecip, col = NDVI), size = 1)+
scale_fill_viridis_c(aesthetics = c("fill", "col"), name = "NDVI")+
labs(title = "Contour plot SELM")+
theme_minimal(base_size = 10)
mtrx3d <- predict(poly.ndvi, newdata = data.fit, type = "response")
mtrx.melt <- cbind(data.fit, mtrx3d)
names(mtrx.melt) <- c("Temp","logPrecip","NDVI")
cp2 <- ggplot(mtrx.melt,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = fennoscandia, aes(x=Temp, y = logPrecip, col = NDVI), size = 1)+
scale_fill_viridis_c(aesthetics = c("fill", "col"), name = "NDVI")+
labs(title = "Contour plot polynomial model")+
theme_minimal(base_size = 10)
mtrx3d <- predict(glm.ndvi, newdata = data.fit, type = "response")
mtrx.melt <- cbind(data.fit, mtrx3d)
names(mtrx.melt) <- c("Temp","logPrecip","NDVI")
cp3 <- ggplot(mtrx.melt,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = fennoscandia, aes(x=Temp, y = logPrecip, col = NDVI), size = 1)+
scale_fill_viridis_c(aesthetics = c("fill", "col"), name = "NDVI")+
labs(title = "Contour plot binomial model")+
theme_minimal(base_size = 10)
mtrx3d <- predict(betareg.ndvi, newdata = data.fit, type = "response")
mtrx.melt <- cbind(data.fit, mtrx3d)
names(mtrx.melt) <- c("Temp","logPrecip","NDVI")
cp4 <- ggplot(mtrx.melt,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..), col = "black")+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = fennoscandia, aes(x=Temp, y = logPrecip, col = NDVI), size = 1)+
scale_fill_viridis_c(aesthetics = c("fill", "col"), name = "NDVI")+
labs(title = "Contour plot beta model")+
theme_minimal(base_size = 10)
grid.arrange(cp1,cp2,cp3,cp4, nrow = 2, ncol = 2)
The binomial model preforms better: low residuals and prediction range between 0 and 1. This model will be used to forecast NDVI in 2050 and in 2100.
df.ssp126 <- readRDS("WR/df.ssp126.rds")
df.ssp370 <- readRDS("WR/df.ssp370.rds")
df.ssp126.2100 <- readRDS("WR/df.ssp126.2100.rds")
df.ssp370.2100 <- readRDS("WR/df.ssp370.2100.rds")
ndvi.model <- readRDS("WR/betareg.ndvi.rds")
wr.kmat <- readRDS("WR/wr.kmat.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- wr.sf.95$Precip %>% log()
#NDVI 126 2050
# ndvi.ssp126.fit <- predict(object = ndvi.model, newdata = df.ssp126, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% as.data.frame() %>% setNames(c("ndvi","trend","signal"))
# ndvi.ssp126 <- cbind(df.ssp126, ndvi.ssp126.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
ndvi.ssp126.fit <- predict(object = ndvi.model, newdata = df.ssp126, type = "response") %>% as.data.frame() %>% setNames("NDVI")
ndvi.ssp126 <- cbind(df.ssp126, ndvi.ssp126.fit)
saveRDS(ndvi.ssp126,"WR/ndvi.ssp126.rds")
# NDVI 370 2050
# ndvi.ssp370.fit <- predict(object = ndvi.model, newdata = df.ssp370, listw = wr.kmat, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% as.data.frame() %>% setNames(c("ndvi","trend","signal"))
#
# ndvi.ssp370 <- cbind(df.ssp370, ndvi.ssp370.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
ndvi.ssp370.fit <- predict(object = ndvi.model, newdata = df.ssp370, type = "response") %>% as.data.frame() %>% setNames("NDVI")
ndvi.ssp370 <- cbind(df.ssp370, ndvi.ssp370.fit)
saveRDS(ndvi.ssp370,"WR/ndvi.ssp370.rds")
#NDVI 126 2010
# ndvi.ssp126.2100.fit <- predict(object = ndvi.model, newdata = df.ssp126.2100, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% as.data.frame() %>% setNames(c("ndvi","trend","signal"))
#
# ndvi.ssp126.2100 <- cbind(df.ssp126.2100, ndvi.ssp126.2100.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
ndvi.ssp126.2100.fit <- predict(object = ndvi.model, newdata = df.ssp126.2100, type = "response") %>% as.data.frame() %>% setNames("NDVI")
ndvi.ssp126.2100 <- cbind(df.ssp126.2100, ndvi.ssp126.2100.fit)
saveRDS(ndvi.ssp126.2100,"WR/ndvi.ssp126.2100.rds")
#NDVI 370 2100
# ndvi.ssp370.2100.fit <- predict(object = ndvi.model, newdata = df.ssp370.2100, listw = wr.kmat, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% setNames(c("ndvi","trend","signal"))
# ndvi.ssp370.2100 <- cbind(df.ssp370.2100, ndvi.ssp370.2100.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
ndvi.ssp370.2100.fit <- predict(object = ndvi.model, newdata = df.ssp370.2100, type = "response") %>% as.data.frame() %>% setNames("NDVI")
ndvi.ssp370.2100 <- cbind(df.ssp370.2100, ndvi.ssp370.2100.fit)
saveRDS(ndvi.ssp370.2100,"WR/ndvi.ssp370.2100.rds")
In the historical dataset in Worldclim, the runoff is given as a yearly average in kg/m2/s, or mm/s.
In the CNRM-CM6 forecast we get the monthly forecast for runoff (in mm/s) and convert it in mm/y.
listgrob <- list("CORDEX_future/runoff.ssp126.2050.png", "CORDEX_future/runoff.ssp370.2050.png","CORDEX_future/runoff.ssp126.2100.png","CORDEX_future/runoff.ssp370.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
ggarrange(plotlist = listgrob, ncol = 2, nrow = 2, labels = c("SSP126-2050","SSP370-2050","SSP126-2100","SSP370-2100"))
wr.sf <- readRDS("WR/wr.sf.rds")
ssp126.runoff.wr <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp126_r1i1p1f2_gr_201501-210012.nc", bands = c(313:552))
ssp126.runoff.wr.mean <- ssp126.runoff.wr %>% mean()
ssp126.runoff.wr.mean.crop <- raster::crop(ssp126.runoff.wr.mean,c(0,35,55,73))
plot(ssp126.runoff.wr.mean.crop)
ssp126.runoff <- raster::extract(ssp126.runoff.wr.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp126.runoff) <- c("gid","runoff")
ssp126.runoff$Runoff <- ssp126.runoff$runoff *365*24*60*60
saveRDS(ssp126.runoff,"WR/ssp126.runoff.2050.rds")
ssp370.runoff.wr <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp370_r1i1p1f2_gr_201501-210012.nc", bands = c(313:552))
ssp370.runoff.wr.mean <- ssp370.runoff.wr %>% mean()
ssp370.runoff.wr.mean.crop <- raster::crop(ssp370.runoff.wr.mean,c(0,35,55,73))
plot(ssp370.runoff.wr.mean.crop)
ssp370.runoff <- raster::extract(ssp370.runoff.wr.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp370.runoff) <- c("gid","runoff")
ssp370.runoff$Runoff <- ssp370.runoff$runoff *365*24*60*60
saveRDS(ssp370.runoff,"WR/ssp370.runoff.2050.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
raster.ssp126.runoff.2100 <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp126_r1i1p1f2_gr_201501-210012.nc", bands = c(805:1032))
raster.ssp126.runoff.2100.mean <- raster.ssp126.runoff.2100 %>% mean()
raster.ssp126.runoff.2100.mean.crop <- raster::crop(raster.ssp126.runoff.2100.mean,c(0,35,55,73))
plot(raster.ssp126.runoff.2100.mean.crop)
ssp126.runoff.2100 <- raster::extract(raster.ssp126.runoff.2100.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp126.runoff.2100) <- c("gid","runoff")
ssp126.runoff.2100$Runoff <- ssp126.runoff.2100$runoff *365*24*60*60
saveRDS(ssp126.runoff.2100,"WR/ssp126.runoff.2100.rds")
raster.ssp370.runoff.2100 <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp370_r1i1p1f2_gr_201501-210012.nc", bands = c(805:1032))
raster.ssp370.runoff.2100.mean <- raster.ssp370.runoff.2100 %>% mean()
raster.ssp370.runoff.2100.mean.crop <- raster::crop(raster.ssp370.runoff.2100.mean,c(0,35,55,73))
plot(raster.ssp370.runoff.2100.mean.crop)
ssp370.runoff.2100 <- raster::extract(raster.ssp370.runoff.2100.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp370.runoff.2100) <- c("gid","runoff")
ssp370.runoff.2100$Runoff <- ssp370.runoff.2100$runoff *365*24*60*60
saveRDS(ssp370.runoff.2100,"WR/ssp370.runoff.2100.rds")
Future TNdep will largely depend on the public policies implemented under the scenarios SSP 1-2.6 and 3-7.0.
TNdep data are from 2000.
Emission reduction between 2005 and 2020 (European Environment Agency 2021):
Predictions for future emissions for NOx (summary for policymakers AR6 (Masson-Delmotte et al. 2021)):
Assumption for NH3:
ndep.wr <- readRDS("WR/ndep.wr.rds")
ndep.sum <- ndep.wr@data %>% apply(2, sum)
(ndep.sum["woxn"]+ndep.sum["doxn"])/ndep.sum["tndep"]
(ndep.sum["wrdn"]+ndep.sum["drdn"])/ndep.sum["tndep"]
# SSP126 2050
ssp126.ndep.2050 <- data.frame(gid = ndep.wr$gid)
ssp126.ndep.2050$oxn <- ndep.wr$woxn*0.55*0.8 + ndep.wr$doxn*0.55*0.8 # -45% for decrease between 2000 and 2020, -20% between 2020 and 2050
ssp126.ndep.2050$rdn <- ndep.wr$wrdn*0.9*0.8 + ndep.wr$drdn*0.9*0.8 # -10% for decrease between 2000 and 2020, -20% between 2020 and 2050
ssp126.ndep.2050$TNdep <- ssp126.ndep.2050$oxn + ssp126.ndep.2050$rdn
saveRDS(ssp126.ndep.2050,"WR/ssp126.ndep.2050.rds")
# SSP126 2100
ssp126.ndep.2100 <- data.frame(gid = ndep.wr$gid)
ssp126.ndep.2100$oxn <- ndep.wr$woxn*0.55*0.8+ndep.wr$doxn*0.55*0.8 # -45% for decrease between 2000 and 2020, -20% between 2020 and 2100
ssp126.ndep.2100$rdn <- ndep.wr$wrdn*0.9*0.8+ndep.wr$drdn*0.9*0.8 # -10% for decrease between 2000 and 2020, -20% between 2020 and 2100
ssp126.ndep.2100$TNdep <- ssp126.ndep.2100$oxn + ssp126.ndep.2100$rdn
saveRDS(ssp126.ndep.2100,"WR/ssp126.ndep.2100.rds")
# SSP370 2050
ssp370.ndep.2050 <- data.frame(gid = ndep.wr$gid)
ssp370.ndep.2050$oxn <- ndep.wr$woxn*0.55*1.37+ndep.wr$doxn*0.55*1.37 # -45% for decrease between 2000 and 2020, 137% increase between 2020 and 2100
ssp370.ndep.2050$rdn <- ndep.wr$wrdn*0.9+ndep.wr$drdn*0.9 # -10% for decrease between 2000 and 2020, no change between 2020 and 2100
ssp370.ndep.2050$TNdep <- ssp370.ndep.2050$oxn + ssp370.ndep.2050$rdn
saveRDS(ssp370.ndep.2050,"WR/ssp370.ndep.2050.rds")
# SSP370 2100
ssp370.ndep.2100 <- data.frame(gid = ndep.wr$gid)
ssp370.ndep.2100$oxn <- ndep.wr$woxn*0.55*1.9+ndep.wr$doxn*0.55*1.9 # -45% for decrease between 2000 and 2020, 190% increase between 2020 and 2100
ssp370.ndep.2100$rdn <- ndep.wr$wrdn*0.9+ndep.wr$drdn*0.9 # -10% for decrease between 2000 and 2020, no change between 2020 and 2100
ssp370.ndep.2100$TNdep <- ssp370.ndep.2100$oxn + ssp370.ndep.2100$rdn
saveRDS(ssp370.ndep.2100,"WR/ssp370.ndep.2100.rds")
The datasets were gathered into a single dataframe.
ndvi.ssp126 <- readRDS("WR/ndvi.ssp126.rds")
ndvi.ssp370 <- readRDS("WR/ndvi.ssp370.rds")
ssp126.runoff.2050 <- readRDS("WR/ssp126.runoff.2050.rds")
ssp370.runoff.2050 <- readRDS("WR/ssp370.runoff.2050.rds")
ssp126.ndep.2050 <- readRDS("WR/ssp126.ndep.2050.rds")
ssp370.ndep.2050 <- readRDS("WR/ssp370.ndep.2050.rds")
ndvi.ssp126.2100 <- readRDS("WR/ndvi.ssp126.2100.rds")
ndvi.ssp370.2100 <- readRDS("WR/ndvi.ssp370.2100.rds")
ssp126.runoff.2100 <- readRDS("WR/ssp126.runoff.2100.rds")
ssp370.runoff.2100 <- readRDS("WR/ssp370.runoff.2100.rds")
ssp126.ndep.2100 <- readRDS("WR/ssp126.ndep.2100.rds")
ssp370.ndep.2100 <- readRDS("WR/ssp370.ndep.2100.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
# 2050 -----
wr.sf.126.2050 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(select(ndvi.ssp126, !c("Arable", "Runoff")), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp126.runoff.2050)), by = "gid") %>% merge(ssp126.ndep.2050)
wr.sf.370.2050 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(select(ndvi.ssp370, !c("Arable","Runoff")), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp370.runoff.2050)), by = "gid") %>% merge(ssp370.ndep.2050)
na.lines.126 <- which(complete.cases(st_drop_geometry(wr.sf.126.2050))==F)
na.lines.370 <- which(complete.cases(st_drop_geometry(wr.sf.370.2050))==F)
wr.sf.126.2050 <- wr.sf.126.2050[-na.lines.126,]
wr.sf.370.2050 <- wr.sf.370.2050[-na.lines.370,]
wr.sf.126.2050$logRunoff <- wr.sf.126.2050$Runoff %>% log()
wr.sf.370.2050$logRunoff <- wr.sf.370.2050$Runoff %>% log()
saveRDS(wr.sf.126.2050, "WR/wr.sf.126.2050.rds")
saveRDS(wr.sf.370.2050, "WR/wr.sf.370.2050.rds")
# 2100 ----
wr.sf.126.2100 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(select(ndvi.ssp126.2100, !c("Arable","Runoff")), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp126.runoff.2100)), by = "gid") %>% merge(ssp126.ndep.2100)
wr.sf.370.2100 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(select(ndvi.ssp370.2100, !c("Arable","Runoff")), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp370.runoff.2100)), by = "gid") %>% merge(ssp370.ndep.2100)
na.lines.126 <- which(complete.cases(st_drop_geometry(wr.sf.126.2100))==F)
na.lines.370 <- which(complete.cases(st_drop_geometry(wr.sf.370.2100))==F)
wr.sf.126.2100 <- wr.sf.126.2100[-na.lines.126,]
wr.sf.370.2100 <- wr.sf.370.2100[-na.lines.370,]
wr.sf.126.2100$logRunoff <- wr.sf.126.2100$Runoff %>% log()
wr.sf.370.2100$logRunoff <- wr.sf.370.2100$Runoff %>% log()
saveRDS(wr.sf.126.2100, "WR/wr.sf.126.2100.rds")
saveRDS(wr.sf.370.2100, "WR/wr.sf.370.2100.rds")
Then we plot the difference between the forecasted variables in 2050 and 2100, compared to 1995.
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
t1 <- ggplot(wr.sf.126.2050)+geom_point(aes(x = Temp95, y = Temp))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Temperature SSP 1-2.6, 2050", x = "Temperature 1995 (°C)")+theme_minimal()
t2 <- ggplot(wr.sf.126.2100)+geom_point(aes(x = Temp95, y = Temp))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Temperature SSP 1-2.6, 2100", x = "Temperature 1995 (°C)")+theme_minimal()
t3 <- ggplot(wr.sf.370.2050)+geom_point(aes(x = Temp95, y = Temp))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Temperature SSP 3-7.0, 2050", x = "Temperature 1995 (°C)")+theme_minimal()
t4 <- ggplot(wr.sf.370.2100)+geom_point(aes(x = Temp95, y = Temp))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Temperature SSP 3-7.0, 2100", x = "Temperature 1995 (°C)")+theme_minimal()
ggarrange(t1,t2,t3,t4, nrow = 2, ncol = 2)
p1 <- ggplot(wr.sf.126.2050)+geom_point(aes(x = Prec95, y = Precip))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Precipitations SSP 1-2.6, 2050", x = "Precipitations 1995 (mm/y)")+theme_minimal()
p2 <- ggplot(wr.sf.126.2100)+geom_point(aes(x = Prec95, y = Precip))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y= "Precipitations SSP 1-2.6, 2100", x = "Precipitations 1995 (mm/y)")+theme_minimal()
p3 <- ggplot(wr.sf.370.2050)+geom_point(aes(x = Prec95, y = Precip))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Precipitations SSP 3-7.0, 2050", x = "Precipitations 1995 (mm/y)")+theme_minimal()
p4 <- ggplot(wr.sf.370.2100)+geom_point(aes(x = Prec95, y = Precip))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Precipitations SSP 3-7.0, 2100", x = "Precipitations 1995 (mm/y)")+theme_minimal()
ggarrange(p1,p2,p3,p4, nrow = 2, ncol = 2)
# d1 <- ggplot(wr.sf.126.2050)+geom_point(aes(x = TNdep95, y = TNdep))+geom_abline(slope = 1, intercept = 0, col = "red")+theme_minimal()
# d2 <- ggplot(wr.sf.126.2100)+geom_point(aes(x = TNdep95, y = TNdep))+geom_abline(slope = 1, intercept = 0, col = "red")+theme_minimal()
# d3 <- ggplot(wr.sf.370.2050)+geom_point(aes(x = TNdep95, y = TNdep))+geom_abline(slope = 1, intercept = 0, col = "red")+theme_minimal()
# d4 <- ggplot(wr.sf.370.2100)+geom_point(aes(x = TNdep95, y = TNdep))+geom_abline(slope = 1, intercept = 0, col = "red")+theme_minimal()
#
# ggarrange(d1,d2,d3,d4, nrow = 2, ncol = 2)
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
n1 <- ggplot(wr.sf.126.2050)+geom_point(aes(x = NDVI95, y = NDVI))+ geom_abline(slope = 1, intercept = 0, col = "red")+
ylim(0,1)+labs(y = "NDVI SSP 1-2.6, 2050")+ theme_minimal()
n2 <- ggplot(wr.sf.126.2100)+geom_point(aes(x = NDVI95, y = NDVI))+geom_abline(slope = 1, intercept = 0, col = "red")+
ylim(0,1)+labs(y = "NDVI SSP 1-2.6, 2100")+theme_minimal()
n3 <- ggplot(wr.sf.370.2050)+geom_point(aes(x = NDVI95, y = NDVI))+geom_abline(slope = 1, intercept = 0, col = "red")+
ylim(0,1)+labs(y = "NDVI SSP 3-7.0, 2050")+theme_minimal()
n4 <- ggplot(wr.sf.370.2100)+geom_point(aes(x = NDVI95, y = NDVI))+geom_abline(slope = 1, intercept = 0, col = "red")+
ylim(0,1)+labs(y = "NDVI SSP 3-7.0, 2100")+theme_minimal()
ggarrange(n1,n2,n3,n4, nrow = 2, ncol = 2)
glm.ndvi <- readRDS("WR/glm.ndvi.rds")
fennoscandia <- readRDS("fennoscandia.33.rds")
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
par(mfrow = c(2,2))
xtemp1 <- seq(min(wr.sf.126.2050$Temp), max(wr.sf.126.2050$Temp), 0.1)
yprecip1 <- seq(min(wr.sf.126.2050$logPrecip), max(wr.sf.126.2050$logPrecip), 0.1)
data.fit1 <- expand.grid(Temp = xtemp1, logPrecip = yprecip1)
mtrx3d1 <- predict(glm.ndvi, newdata = data.fit1, type = "response")
mtrx.melt1 <- cbind(data.fit1, mtrx3d1)
names(mtrx.melt1) <- c("Temp","logPrecip","NDVI")
fcp1 <- ggplot(mtrx.melt1,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = wr.sf.126.2050, aes(x=Temp, y = logPrecip), size = 1)+
scale_fill_viridis_c()+
labs(title = "Contour plot NDVI SSP 1-2.6, 2050")
xtemp2 <- seq(min(wr.sf.126.2100$Temp), max(wr.sf.126.2100$Temp), 0.1)
yprecip2 <- seq(min(wr.sf.126.2100$logPrecip), max(wr.sf.126.2100$logPrecip), 0.1)
data.fit2 <- expand.grid(Temp = xtemp2, logPrecip = yprecip2)
mtrx3d2 <- predict(glm.ndvi, newdata = data.fit2, type = "response")
mtrx.melt2 <- cbind(data.fit2, mtrx3d2)
names(mtrx.melt2) <- c("Temp","logPrecip","NDVI")
fcp2 <- ggplot(mtrx.melt2,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = wr.sf.126.2100, aes(x=Temp, y = logPrecip), size = 1)+
scale_fill_viridis_c()+
labs(title = "Contour plot NDVI SSP 1-2.6, 2100")
xtemp3 <- seq(min(wr.sf.370.2050$Temp), max(wr.sf.370.2050$Temp), 0.1)
yprecip3 <- seq(min(wr.sf.370.2050$logPrecip), max(wr.sf.370.2050$logPrecip), 0.1)
data.fit3 <- expand.grid(Temp = xtemp3, logPrecip = yprecip3)
mtrx3d3 <- predict(glm.ndvi, newdata = data.fit3, type = "response")
mtrx.melt3 <- cbind(data.fit3, mtrx3d3)
names(mtrx.melt3) <- c("Temp","logPrecip","NDVI")
fcp3 <- ggplot(mtrx.melt3,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = wr.sf.370.2050, aes(x=Temp, y = logPrecip), size = 1)+
scale_fill_viridis_c()+
labs(title = "Contour plot NDVI SSP 3-7.0, 2050")
xtemp4 <- seq(min(wr.sf.370.2100$Temp), max(wr.sf.370.2100$Temp), 0.1)
yprecip4 <- seq(min(wr.sf.370.2100$logPrecip), max(wr.sf.370.2100$logPrecip), 0.1)
data.fit4 <- expand.grid(Temp = xtemp4, logPrecip = yprecip4)
mtrx3d4 <- predict(glm.ndvi, newdata = data.fit4, type = "response")
mtrx.melt4 <- cbind(data.fit4, mtrx3d4)
names(mtrx.melt4) <- c("Temp","logPrecip","NDVI")
fcp4 <- ggplot(mtrx.melt4,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = wr.sf.370.2100, aes(x=Temp, y = logPrecip), size = 1)+
scale_fill_viridis_c()+
labs(title = "Contour plot NDVI SSP 3-7.0, 2100")
grid.arrange(fcp1,fcp2,fcp3,fcp4,nrow=2,ncol=2)
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
mn0 <- ggplot(wr.sf.126.2050)+geom_sf(aes(fill = NDVI95, col = NDVI95)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI 1995", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn0, filename = "WR/NDVI_1995.png",dpi = "retina", width = 12, height = 16, units = "cm")
mn1 <- ggplot(wr.sf.126.2050)+geom_sf(aes(fill = NDVI, col = NDVI)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI SSP126-2050", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn1, filename = "WR/NDVI_SSP126_2050.png",dpi = "retina", width = 12, height = 16, units = "cm")
mn2 <- ggplot(wr.sf.126.2100)+geom_sf(aes(fill = NDVI, col = NDVI)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI SSP126-2100", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn2, filename = "WR/NDVI_SSP126_2100.png",dpi = "retina", width = 12, height = 16, units = "cm")
mn3 <- ggplot(wr.sf.370.2050)+geom_sf(aes(fill = NDVI, col = NDVI)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI SSP370-2050", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn3, filename = "WR/NDVI_SSP370_2050.png",dpi = "retina", width = 12, height = 16, units = "cm")
mn4 <- ggplot(wr.sf.370.2100)+geom_sf(aes(fill = NDVI, col = NDVI)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI SSP370-2100", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn4, filename = "WR/NDVI_SSP370_2100.png",dpi = "retina", width = 12, height = 16, units = "cm")
listgrob <- list("WR/NDVI_1995.png","white.png","WR/NDVI_SSP126_2050.png", "WR/NDVI_SSP126_2100.png","WR/NDVI_SSP370_2050.png","WR/NDVI_SSP370_2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 3, align = "v", greedy = T)
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
diff.var <- function(sf, vars1, vars2){
sf2 <- sf
for (i in 1:length(vars2)){
newname <- paste("diff", vars2[i], sep = "")
x <- sf2[,vars2[i]] %>% st_drop_geometry()
y <- sf2[,vars1[i]] %>% st_drop_geometry()
z <- x-y
names(z) <- newname
sf2 <- cbind(sf2, z)
}
return(sf2)
}
wr.sf.126.2050 <- diff.var(wr.sf.126.2050, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
wr.sf.370.2050 <- diff.var(wr.sf.370.2050, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
wr.sf.126.2100 <- diff.var(wr.sf.126.2100, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
wr.sf.370.2100 <- diff.var(wr.sf.370.2100, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
list.sf <- list(wr.sf.126.2050, wr.sf.126.2100, wr.sf.370.2050, wr.sf.370.2100)
names(list.sf) <- c("SSP126.2050", "SSP126.2100", "SSP370.2050", "SSP370.2100")
all.diff <- data.frame(gid = c(wr.sf.126.2050$gid, wr.sf.370.2050$gid, wr.sf.126.2100$gid, wr.sf.370.2100$gid))
for(x in c("diffTemp","diffPrecip","diffNDVI","diffRunoff","diffTNdep")){
all.vec <- c()
for(i in 1:length(list.sf)){
vec <- select(st_drop_geometry(list.sf[[i]]),all_of(x))
all.vec <<- rbind(all.vec,vec)
}
all.diff <- cbind(all.diff, all.vec)
}
max.lim <- apply(all.diff, 2, max, na.rm = T)
min.lim <- apply(all.diff, 2, min, na.rm = T)
# x1 <- filter(wr.sf.126.2050, NDVI95 > 0.4)
# x2 <- filter(wr.sf.370.2050, NDVI95 > 0.4)
# x3 <- filter(wr.sf.126.2100, NDVI95 > 0.4)
# x4 <- filter(wr.sf.370.2100, NDVI95 > 0.4)
# all.ndvi <- c(x1$NDVI/x1$NDVI95*100-100,x2$NDVI/x2$NDVI95*100-100, x3$NDVI/x3$NDVI95*100-100, x4$NDVI/x4$NDVI95*100-100)
# palettes <- data.frame(diffTemp = brewer.pal(name = "RdYlBu", n = 9)[c(1,5,9)],
# diffPrecip = brewer.pal(name = "RdYlBu", n = 9)[c(9,5,1)],
# diffNDVI = brewer.pal(name = "RdYlGn", n = 9)[c(9,5,1)],
# diffRunoff = brewer.pal(name = "RdYlBu", n = 9)[c(9,5,1)],
# diffTNdep = brewer.pal(name = "RdYlBu", n = 9)[c(9,5,1)])
legends <- c("Difference Temperature (°C)", "Difference Precipitation (mm/y)", "Difference NDVI", "Difference Runoff (mm/y)","DIfference Nitrogen Deposition (ug/m2)")
diffs <- c("diffTemp","diffPrecip","diffNDVI","diffRunoff","diffTNdep")
rev.pal <- c(T,F,T,F,T)
for(i in c(1:length(list.sf))){
for(x in diffs[3]){
g <- ggplot(list.sf[[i]])+geom_sf(aes(fill = .data[[x]], col = .data[[x]])) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = rev.pal[grep(x, diffs)] , aesthetics = c("fill","col"), name = paste(legends[grep(x, diffs)],"\n",names(list.sf)[i], sep = ""), limits = c(min.lim[[x]], max.lim[[x]])) +
# scale_fill_gradient2(low=palettes[3, grep(x,names(palettes))],mid=palettes[2, grep(x,names(palettes))],high = palettes[1, grep(x,names(palettes))],
# midpoint = 0, aesthetics = c("fill","col"),
# name = paste(legends[grep(x, diffs)],"\n",names(list.sf)[i], sep = ""), limits = c(min.lim[[x]], max.lim[[x]])) +
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = g, filename = paste("WR/", names(list.sf)[i], ".", x, ".png",sep = ""),dpi = "retina", width = 12, height = 16, units = "cm")
}
}
white <- ggplot()+theme_void()
ggsave(plot = white, filename = "white.png",dpi = "retina", width = 12, height = 16, units = "cm")
listgrob <- list("WR/SSP126.2050.diffTemp.png", "WR/SSP126.2100.diffTemp.png","WR/SSP370.2050.diffTemp.png","WR/SSP370.2100.diffTemp.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
listgrob <- list("WR/SSP126.2050.diffPrecip.png","WR/SSP126.2100.diffPrecip.png","WR/SSP370.2050.diffPrecip.png","WR/SSP370.2100.diffPrecip.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
listgrob <- list("WR/SSP126.2050.diffNDVI.png","WR/SSP126.2100.diffNDVI.png","WR/SSP370.2050.diffNDVI.png","WR/SSP370.2100.diffNDVI.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
listgrob <- list("WR/SSP126.2050.diffRunoff.png","WR/SSP126.2100.diffRunoff.png","WR/SSP370.2050.diffRunoff.png","WR/SSP370.2100.diffRunoff.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
listgrob <- list("WR/SSP126.2050.diffTNdep.png","WR/SSP126.2100.diffTNdep.png","WR/SSP370.2050.diffTNdep.png","WR/SSP370.2100.diffTNdep.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
The average concentration of total organic carbon in the coastal drainage basins was predicted based on the SELM with 5 predictors, fitted on the Fennoscandian dataset. The difference between the predicted log(TOC) for 2050 and 2100 under the SSP 1-2.6 and 3-7.0 scenarios, compared to the log(TOC) predicted for 1995, was computed and plotted.
fennoscandia <- readRDS("fennoscandia.33.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
wr.ssp126.2050 <- readRDS("WR/wr.ssp126.2050.rds")
wr.ssp370.2050 <- readRDS("WR/wr.ssp370.2050.rds")
wr.ssp126.2100 <- readRDS("WR/wr.ssp126.2050.rds")
wr.ssp370.2100 <- readRDS("WR/wr.ssp370.2050.rds")
temperature <- c(fennoscandia$Temp, sem.pred.wr$Temp,wr.ssp126.2050$Temp,wr.ssp370.2050$Temp, wr.ssp126.2100$Temp,wr.ssp370.2100$Temp)
precipitation <- c(fennoscandia$Precip, sem.pred.wr$Precip,wr.ssp126.2050$Precip,wr.ssp370.2050$Precip, wr.ssp126.2100$Precip,wr.ssp370.2100$Precip)
ndvi <- c(fennoscandia$NDVI, sem.pred.wr$NDVI,wr.ssp126.2050$NDVI,wr.ssp370.2050$NDVI, wr.ssp126.2100$NDVI,wr.ssp370.2100$NDVI)
runoff <- c(fennoscandia$Runoff, sem.pred.wr$Runoff,wr.ssp126.2050$Runoff,wr.ssp370.2050$Runoff, wr.ssp126.2100$Runoff,wr.ssp370.2100$Runoff)
tndep <- c(fennoscandia$TNdep, sem.pred.wr$TNdep,wr.ssp126.2050$TNdep,wr.ssp370.2050$TNdep, wr.ssp126.2100$TNdep,wr.ssp370.2100$TNdep)
toc <- c(fennoscandia$TOC, sem.pred.wr$TOC95,wr.ssp126.2050$TOC,wr.ssp370.2050$TOC, wr.ssp126.2100$TOC,wr.ssp370.2100$TOC)
neworder <- c("Fennoscandia","WR1995","SSP126.2050","SSP126.2100","SSP370.2050","SSP370.2100")
model <- factor(c(rep("Fennoscandia",dim(fennoscandia)[1]),
rep("WR1995",dim(sem.pred.wr)[1]),
rep("SSP126.2050", dim(wr.ssp126.2050)[1]),
rep("SSP370.2050", dim(wr.ssp370.2050)[1]),
rep("SSP126.2100", dim(wr.ssp126.2100)[1]),
rep("SSP370.2100", dim(wr.ssp370.2100)[1])),
levels = neworder)
compare.ssp <- data.frame(value = c(temperature, precipitation, runoff, ndvi, toc), variable = c(rep("Temperature (°C)", length(temperature)), rep("Precipitation (mm)", length(precipitation)), rep("Runoff (mm/y)",length(runoff)), rep("NDVI",length(ndvi)), rep("TOC (mg/L)", length(toc))), model = rep(model,5))
ggplot(compare.ssp)+geom_boxplot(aes(y= value,col=variable))+
labs(y="", x = "")+scale_color_viridis_d(end = 0.8)+
theme_minimal()+
theme(legend.position = "", axis.text.x = element_blank())+
facet_grid(cols = vars(model), rows = vars(factor(variable, levels = c("Temperature (°C)","Precipitation (mm)","Runoff (mm/y)","NDVI","TOC (mg/L)"))), scales = "free_y", shrink = T)
ggsave("WR/compare.ssp.png", width = 30, height = 25, units = "cm", dpi = "retina")
knitr::include_graphics("WR/compare.ssp.png")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp126 <- predict(sem.fennoscandia.5,wr.sf.126.2050,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp126 <- cbind(wr.sf.126.2050,wr.pred.ssp126)
wr.ssp126$TOC <- exp(wr.ssp126$fit)
wr.ssp126 <- merge(wr.ssp126, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp126,"WR/wr.ssp126.2050.rds")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp370 <- predict(sem.fennoscandia.5,wr.sf.370.2050,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp370 <- cbind(wr.sf.370.2050,wr.pred.ssp370)
wr.ssp370$TOC <- exp(wr.ssp370$fit)
wr.ssp370 <- merge(wr.ssp370, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp370,"WR/wr.ssp370.2050.rds")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp126 <- predict(sem.fennoscandia.5,wr.sf.126.2100,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp126 <- cbind(wr.sf.126.2100,wr.pred.ssp126)
wr.ssp126$TOC <- exp(wr.ssp126$fit)
wr.ssp126 <- merge(wr.ssp126, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp126,"WR/wr.ssp126.2100.rds")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp370 <- predict(sem.fennoscandia.5,wr.sf.370.2100,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp370 <- cbind(wr.sf.370.2100,wr.pred.ssp370)
wr.ssp370$TOC <- exp(wr.ssp370$fit)
wr.ssp370 <- merge(wr.ssp370, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp370,"WR/wr.ssp370.2100.rds")
The average TOC concentration was plotted.
wr.ssp126.2050 <- readRDS("WR/wr.ssp126.2050.rds")
wr.ssp370.2050 <- readRDS("WR/wr.ssp370.2050.rds")
wr.ssp126.2100 <- readRDS("WR/wr.ssp126.2100.rds")
wr.ssp370.2100 <- readRDS("WR/wr.ssp370.2100.rds")
list.sf <- list(wr.ssp126.2050, wr.ssp126.2100, wr.ssp370.2050, wr.ssp370.2100)
names(list.sf) <- c("SSP126.2050", "SSP126.2100", "SSP370.2050", "SSP370.2100")
for(i in names(list.sf)){
list.sf[[i]]$diffTOC <- list.sf[[i]]$TOC - list.sf[[i]]$TOC95
list.sf[[i]]$diff.log <- list.sf[[i]]$fit - list.sf[[i]]$wr.pred.fit
list.sf[[i]]$exp.diff.log <- exp(list.sf[[i]]$diff.log)
list.sf[[i]]$percent.change <- (list.sf[[i]]$exp.diff.log - 1) * 100
}
all.diff <- data.frame(gid = c(wr.ssp126.2050$gid, wr.ssp370.2050$gid, wr.ssp126.2100$gid, wr.ssp370.2100$gid))
for(x in c("fit","diffTOC","diff.log","exp.diff.log","percent.change")){
all.vec <- c()
for(i in names(list.sf)){
vec <- dplyr::select(st_drop_geometry(list.sf[[i]]),all_of(x))
all.vec <<- rbind(all.vec,vec)
}
all.diff <- cbind(all.diff, all.vec)
}
max.lim <- apply(all.diff, 2, max, na.rm = T)
min.lim <- apply(all.diff, 2, min, na.rm = T)
t0 <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=wr.pred.fit, col= wr.pred.fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "logTOC in 1995")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave("WR/toc_1995.png", t0 ,dpi = "retina", width = 12, height = 16, units = "cm")
t1 <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=fit, col= fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "a) Forecast logTOC\n SSP1-2.6 - 2050")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
# scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="a) Forecast logTOC\n SSP1-2.6 - 2050", mid = (max.lim[2]-min.lim[2])/2, limits = c(min.lim[2],max.lim[2]), rev = T) +
# theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave("WR/logTOC_ssp126.2050.png", t1 ,dpi = "retina", width = 12, height = 16, units = "cm")
t2 <- ggplot(list.sf[["SSP126.2100"]])+
geom_sf(aes(fill=fit, col= fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "b) Forecast logTOC\n SSP1-2.6 - 2100")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave("WR/logTOC_ssp126.2100.png", t2 ,dpi = "retina", width = 12, height = 16, units = "cm")
t3 <- ggplot(list.sf[["SSP370.2050"]])+
geom_sf(aes(fill=fit, col= fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "c) Forecast logTOC\n SSP3-7.0 - 2050")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave("WR/logTOC_ssp370.2050.png", t3 ,dpi = "retina", width = 12, height = 16, units = "cm")
t4 <- ggplot(list.sf[["SSP370.2100"]])+
geom_sf(aes(fill=fit, col= fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "d) Forecast logTOC\n SSP3-7.0 - 2100")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave("WR/logTOC_ssp370.2100.png", t1 ,dpi = "retina", width = 12, height = 16, units = "cm")
knitr::include_graphics("WR/toc_1995.png")
list.grob <- c("WR/logTOC_ssp126.2050.png", "WR/logTOC_ssp126.2100.png","WR/logTOC_ssp370.2050.png","WR/logTOC_ssp370.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow = 2, ncol = 2, labels = c("a)","b)","c)","d)"), label_size = 10)
p1 <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
# scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],max.lim[5]),
# name = "a) Forecasted percent increase in [TOC]\n SSP1-2.6 - 2050")+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="a) Forecast [TOC] increase in %\n SSP1-2.6 - 2050", mid = 1, limits = c(min.lim[5],max.lim[5]), rev = T) +
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
p2 <- ggplot(list.sf[["SSP126.2100"]])+
geom_sf(aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="b) Forecast [TOC] increase in %\n SSP1-2.6 - 2100", mid = 1, limits = c(min.lim[5],max.lim[5]), rev = T) +
# scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],max.lim[5]),
# name = "b) Forecasted percent change TOC\n SSP1-2.6 - 2100")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
p3 <- ggplot(list.sf[["SSP370.2050"]])+
geom_sf(aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="c) Forecast [TOC] increase in %\n SSP3-7.0 - 2050", mid = 1, limits = c(min.lim[5],max.lim[5]), rev = T) +
# scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],max.lim[5]),
# name = "c) Forecasted percent change TOC\n SSP3-7.0 - 2050")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
p4 <- ggplot(list.sf[["SSP370.2100"]])+
geom_sf(aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="d) Forecast [TOC] increase in %\n SSP3-7.0 - 2100", mid = 1, limits = c(min.lim[5],max.lim[5]), rev = T) +
# scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],max.lim[5]),
# name = "d) Forecasted percent change TOC\n SSP3-7.0 - 2100")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
pp <- grid.arrange(p1, p2, p3, p4, ncol = 2, nrow = 2)
ggsave("forecast.percent.change.toc.png", pp, width = 12, height = 16, units = "cm")
highchangeid <- subset(list.sf$SSP370.2100, exp.diff.log > 25)$gid
pp1 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP126.2050"]], exp.diff.log < 20), aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],20),
name = "a) Forecasted percent change TOC\n SSP1-2.6 - 2050")+
geom_sf(data = subset(list.sf[["SSP126.2050"]], exp.diff.log > 20), aes(fill=exp.diff.log), col= "black", fill = "black", size= 0.05)+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
pp2 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP126.2100"]], exp.diff.log < 20), aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],20),
name = "b) Forecasted percent change TOC\n SSP1-2.6 - 2100")+
geom_sf(data = subset(list.sf[["SSP126.2100"]], exp.diff.log > 20), aes(fill=exp.diff.log), col= "black", fill = "black", size= 0.05)+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
pp3 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP370.2050"]], exp.diff.log < 20), aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],20),
name = "c) Forecasted percent change TOC\n SSP3-7.0 - 2050")+
geom_sf(data = subset(list.sf[["SSP370.2050"]], exp.diff.log > 20), aes(fill=exp.diff.log), col= "black", fill = "black", size= 0.05)+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
pp4 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP370.2100"]], exp.diff.log < 20), aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],20),
name = "d) Forecasted percent change TOC\n SSP3-7.0 - 2100")+
geom_sf(data = subset(list.sf[["SSP370.2100"]], exp.diff.log > 20), aes(fill=exp.diff.log), col= "black", fill = "black", size= 0.05)+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ppp <- grid.arrange(pp1, pp2, pp3, pp4, ncol = 2, nrow = 2)
ggsave("WR/forecast.prop.change.png", ppp, width = 12, height = 16)
c1 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP126.2050"]], percent.change < 307), aes(fill=percent.change, col= percent.change), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="a) Forecast [TOC] change in %\n SSP1-2.6 - 2050", mid = 0, limits = c(min.lim[6],307), rev = T) +
geom_sf(data = subset(list.sf[["SSP126.2050"]], percent.change > 307), fill = "black", col = "black")+
theme_minimal(base_size = 12)+
theme(legend.position = "bottom")
ggsave("WR/percentTOC_ssp126.2050.png", c1 ,dpi = "retina", width = 12, height = 16, units = "cm")
c2 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP126.2100"]], percent.change < 307), aes(fill=percent.change, col= percent.change), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="b) Forecast [TOC] change in %\n SSP1-2.6 - 2100", mid = 0, limits = c(min.lim[6],307), rev = T) +
geom_sf(data = subset(list.sf[["SSP126.2100"]], percent.change > 307), fill = "black", col = "black")+
theme_minimal(base_size = 12)+
theme(legend.position = "bottom")
ggsave("WR/percentTOC_ssp126.2100.png", c2 ,dpi = "retina", width = 12, height = 16, units = "cm")
c3 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP370.2050"]], percent.change < 307), aes(fill=percent.change, col= percent.change), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="c) Forecast [TOC] change in %\n SSP3-7.0 - 2050", mid = 0, limits = c(min.lim[6],307), rev = T) +
geom_sf(data = subset(list.sf[["SSP370.2050"]], percent.change > 307), fill = "black", col = "black")+
theme_minimal(base_size = 12)+
theme(legend.position = "bottom")
ggsave("WR/percentTOC_ssp370.2050.png", c3 ,dpi = "retina", width = 12, height = 16, units = "cm")
c4 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP370.2100"]], percent.change < 307), aes(fill=percent.change, col= percent.change), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="d) Forecast [TOC] change in %\n SSP3-7.0 - 2100", mid = 0, limits = c(min.lim[6],307), rev = T) +
geom_sf(data = subset(list.sf[["SSP370.2100"]], percent.change > 307), fill = "black", col = "black")+
theme_minimal(base_size = 12)+
theme(legend.position = "bottom")
ggsave("WR/percentTOC_ssp370.2100.png", c4 ,dpi = "retina", width = 12, height = 16, units = "cm")
list.grob <- c("WR/percentTOC_ssp126.2050.png", "WR/percentTOC_ssp126.2100.png","WR/percentTOC_ssp370.2050.png","WR/percentTOC_ssp370.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow = 2, ncol = 2, labels = c("a)","b)","c)","d)"), label_size = 10)
g <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="a) Difference logTOC\n SSP1-2.6 - 2050", mid = 0, limits = c(min.lim[4],max.lim[4]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP126.2050.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP126.2100"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="b) Difference logTOC\n SSP1-2.6 - 2100", mid = 0, limits = c(min.lim[4],max.lim[4]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP126.2100.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP370.2050"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="c) Difference logTOC\n SSP3-7.0 - 2050", mid = 0, limits = c(min.lim[4],max.lim[4]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP370.2050.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP370.2100"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="d) Difference logTOC\n SSP3-7.0 - 2100", mid = 0, limits = c(min.lim[4],max.lim[4]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP370.2100.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
listgrob <- list("WR/SSP126.2050.diff.log.png", "WR/SSP126.2100.diff.log.png","WR/SSP370.2050.diff.log.png","WR/SSP370.2100.diff.log.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
ssp <- plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
save_plot("ssp.toc.png",ssp, ncol = 2, nrow = 2, base_height = 4, base_width = 3)
print(ssp)
Most coastal drainage basins in Sweden and Finland drain into to the Baltic Sea (see Supplementary 2), while Norwegian sea watershed drain in Skagerrak, in the North sea, in the Norwegian Sea and in the Barent sea in the North (Sætre 2007). The drainage basin of the Baltic sea matches almost exactly the borders of Sweden and Finland : https://www.grida.no/resources/5324. All contribute to the Norwegian Coastal Current water masses. The estimated exported TOC was computed as:
\[ TOC_exp = [TOC] (mg.L^{-1}) \times Runoff (mm/y, or L.m^{-2}.y^{-1}) \times Watershed area (m^2)\]
Where the TOC concentration (past and future) was obtained by back-transforming the results of the model (log(TOC)) into TOC in mg/L. The result of the above equation, giving the amount of TOC exported in mg, is converted to Tg. Forecast are presented on Figure 19*.
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.centroid <- st_centroid(wr.sf.95, of_largest_polygon = T) %>% dplyr::select("gid","geometry")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
nor <- st_read("Country_shapefile/norway.shp") %>% st_transform(st_crs(wr.centroid))
swe <- st_read("Country_shapefile/sweden.shp") %>% st_transform(st_crs(wr.centroid))
fin <- st_read("Country_shapefile/finland.shp") %>% st_transform(st_crs(wr.centroid))
wr.norge <- wr.centroid[nor,]
wr.norge$Nation <- "Norway"
wr.sweden <- wr.centroid[swe,]
wr.sweden$Nation <- "Sweden"
wr.finland <- wr.centroid[fin,]
wr.finland$Nation <- "Finland"
wr.nation <- rbind(wr.norge,wr.sweden,wr.finland)
wr.sf.95 <- wr.sf.95 %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid") %>% sp::merge(dplyr::select(st_drop_geometry(sem.pred.wr), c("gid","wr.pred.fit","TOC95")), by = "gid")
names(wr.sf.95[which(names(wr.sf.95) == "TOC95")]) <- "TOC"
wr.ssp126.2050 <- readRDS("WR/wr.ssp126.2050.rds")%>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
wr.ssp370.2050 <- readRDS("WR/wr.ssp370.2050.rds") %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
wr.ssp126.2100 <- readRDS("WR/wr.ssp126.2100.rds") %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
wr.ssp370.2100 <- readRDS("WR/wr.ssp370.2100.rds") %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
list.sf <- list(wr.sf.95, wr.ssp126.2050, wr.ssp126.2100, wr.ssp370.2050, wr.ssp370.2100)
names(list.sf) <- c("Data1995","SSP126.2050", "SSP126.2100", "SSP370.2050", "SSP370.2100")
export <- function(sf){
discharge <- sf$Runoff*sf$area # Q in L/y
toc.export <- discharge * sf$TOC * 10^(-15) # amount of exported TOC in Tg/y
sf2 <- cbind(sf,discharge,toc.export)
return(sf2)
}
list.sf <- lapply(list.sf, export)
list.mean <- c()
for(i in 1:length(list.sf)){
list.mean[[i]] <- list.sf[[i]] %>% st_drop_geometry() %>% group_by(Nation) %>% summarise(total = sum(toc.export)) %>% as.data.frame()
}
names(list.mean) <- names(list.sf)
list.diff <- c()
for(i in 1:length(list.mean)){
df1 <- list.mean[[1]]
df2 <- list.mean[[i]]
diff <- df2[,2] - df1[,2]
x <- cbind(list.mean[[i]],diff)
names(x) <- c("Nation", names(list.mean)[i], paste("diff", names(list.mean)[i], sep = ""))
list.diff[[i]] <- x
}
toc.export.df <- list.diff[[1]] %>% merge(list.diff[[2]], by = "Nation") %>% merge(list.diff[[3]], by = "Nation") %>% merge(list.diff[[4]], by = "Nation") %>% merge(list.diff[[5]], by = "Nation") %>% melt(id.vars = c("Nation"), value.name = "TOC_Tg")
toc.export.df$Model <- NA
m.126 <- grep("126",toc.export.df$variable)
m.370 <- grep("370",toc.export.df$variable)
toc.export.df$Model[m.126] <- "SSP126"
toc.export.df$Model[m.370] <- "SSP370"
toc.export.df$Model[-c(m.126,m.370)] <- "Historical"
toc.export.df$Year <- NA
m.95 <- grep("95",toc.export.df$variable)
m.2050 <- grep("2050",toc.export.df$variable)
m.2100 <- grep("2100",toc.export.df$variable)
toc.export.df$Year[m.95] <- "1995"
toc.export.df$Year[m.2050] <- "2050"
toc.export.df$Year[m.2100] <- "2100"
toc.export.df$Diff <- NA
m.d <- grep("diff",toc.export.df$variable)
toc.export.df$Diff[m.d] <- "Difference"
toc.export.df$Diff[-m.d] <- "Absolute"
toc.export.sum <- toc.export.df %>% group_by(Model,Year,Diff) %>% summarise(TOC_Tg = sum(TOC_Tg))
toc.export.sum$Nation <- "Total"
toc.export.df <- rbind(toc.export.sum,toc.export.df)
toc.export.df <- transform(toc.export.df, Nation = factor(Nation, levels = c("Norway","Sweden","Finland","Total")))
toc.export.df$variable <- NULL
toc.export.df$TOC_Tg <- as.numeric(toc.export.df$TOC_Tg)
knitr::kable(toc.export.df) %>% kable_styling()
saveRDS(toc.export.df,"WR/toc.export.df.rds")
ggplot() +
geom_point(data = toc.export.df, aes(x=Year, y = TOC_Tg, col = Model, shape = Model)) +
geom_path(data = toc.export.df, aes(x=Year, y = TOC_Tg, col = Model, group = Model), size = 0.5) +
geom_abline(slope = 0, intercept = 0, col = "firebrick4", size = 0.3)+
scale_color_manual(values = c("firebrick4","steelblue4","darkorange"))+
geom_label(data = filter(toc.export.df, Model == "SSP126" & Diff == "Absolute"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = -0.5, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "SSP126"& Diff == "Difference"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = -0.1, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "SSP370"& Diff == "Absolute"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = 0.5, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "SSP370"& Diff == "Difference"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = 0.1, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "Historical"& Diff == "Absolute"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = 0.5, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "Historical"& Diff == "Difference"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = 0.1, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
labs(x="Year", y= "TOC (Tg)")+
facet_grid(cols=vars(Nation), rows = vars(Diff), scales = "free_y")+
theme_bw(base_size = 12)+
theme(legend.position = "bottom")
ggsave(filename = "WR/TOC_export.png", dpi = "retina", width = 20, height = 15, units = "cm")
toc.export.df <- readRDS("WR/toc.export.df.rds") %>% dplyr::select(c("Nation", "Model", "Year","Diff","TOC_Tg"))
toc.export.df %>% arrange(Model) %>%
knitr::kable(digits = 3,
caption = "TOC export by country",
col.names = c("Country","Model","Year","Absolute/differential export","TOC export (Tg)")) %>%
kable_styling(bootstrap_options = "bordered",
position = "center",
full_width = F)
knitr::include_graphics("WR/TOC_export.png")